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Abstract 

An important class of problems exhibits smooth behaviour on macroscopic space and 
time scales, while only a microscopic evolution law is known. For such time-dependent 
multi-scale problems, an "equation-free" framework has been proposed, of which patch 
dynamics is an essential component. Patch dynamics is designed to perform numerical 
simulations of an unavailable macroscopic equation on macroscopic time and length 
scales; it uses appropriately initialized simulations of the available microscopic model in 
a number of small boxes (patches), which cover only a fraction of the space-time domain. 
To reduce the effect of the artificially introduced box boundaries, we use buffer regions to 
"shield" the boundary artefacts from the interior of the domain for short time intervals. 
We analyze the accuracy of this scheme for a diffusion homogenization problem with 
periodic heterogeneity, and propose a simple heuristic to determine a sufficient buffer 
size. The algorithm performance is illustrated through a set of numerical examples, 
which include a non-linear reaction-diffusion equation and the Kuramoto-Sivashinsky 
equation. 
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1 Introduction 



For an important class of multi-scale problems, a separation of scales prevails between the 
(microscopic, detailed) level of description of the available model, and the (macroscopic, 
continuum) level at which one would like to observe and analyze the system. Consider, 
for example, a kinetic Monte Carlo model of bacterial growth [31]. A stochastic model 
describes the probability of an individual bacterium to run or "tumble", based on the 
rotation of its flagellae. Technically, it would be possible to simply evolve the detailed 
model and observe the macroscopic variables of interest (e.g. cell density), but this could 
be prohibitively expensive. It is known, however, that, under certain conditions, one could 
write a deterministic equation for the evolution of the macroscopic observable (here bacteria 
concentration, the zeroth moment of the evolving distribution) on macroscopic space and 
time scales, but it is hard to obtain an accurate closed formula explicitly. 

The recently proposed equation-free framework [19] can then be used instead of stochas- 
tic time integration in the entire space-time domain. This framework is built around the 
central idea of a coarse time-stepper, which is a time-<5£ map from coarse variables to coarse 
variables. It consists of the following steps: (1) lifting, i.e. the creation of appropriate initial 
conditions for the microscopic model; (2) evolution, using the microscopic model and (pos- 
sibly) some constraints; and (3) restriction, i.e. the projection of the detailed solution to 
the macroscopic observation variables. This coarse time-stepper can subsequently be used 
as "input" for time-stepper based algorithms performing macroscopic numerical analysis 
tasks. These include, for example, time-stepper based bifurcation codes to perform bifur- 
cation analysis for the unavailable macroscopic equation [24, 25, 35, 36]. This approach has 
already been used in several applications [15, 34], and also allows to perform other system 
level tasks, such as control and optimization [33]. 

When dealing with systems that would be described by (in our case, unavailable) par- 
tial differential equations (PDEs), one can also reduce the spatial complexity. For systems 
with one space dimension, the gap-tooth scheme [19] was proposed; it can be generalized in 
several space dimensions. A number of small intervals, separated by large gaps, are intro- 
duced; they qualitatively correspond to mesh points for a traditional, continuum solution of 
the unavailable equation. In higher space dimensions, these intervals would become boxes 
around the coarse mesh points, a term that we will also use throughout this paper. We 
construct a coarse time-<5t map as follows. We first choose a number of macroscopic grid 
points. Then, we choose a small interval around each grid point; initialize the fine scale, 
microscopic solver within each interval consistently with the macroscopic initial condition 
profiles; and provide each box with appropriate boundary conditions. Subsequently, we use 
the microscopic model in each interval to simulate until time 5t, and obtain macroscopic 
information (e.g. by computing the average density in each box) at time 5t. This amounts 
to a coarse time-<5t map; the procedure is then repeated. The resulting scheme has already 
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been used with lattice-Boltzmann simulations of the Fitzhugh-Nagumo dynamics [18, 19] 
and with particle-based simulations of the viscous Burgers equation [10]. 

To increase the efficiency of time integration, one can use the gap-tooth scheme in 
conjuction with any method-of-lines time integration method, such as projective integration 
[8]. We then perform a number of gap-tooth steps of size 5t to obtain an estimate of the 
time derivative of the unavailable macroscopic equation. This estimate is subsequently used 
to perform a time step of size At ^> 5t. This combination has been termed patch dynamics 
[19]. 

In this paper, we will study the patch dynamics scheme for a model diffusion homog- 
enization problem. Here, the microscopic equation is a diffusion equation with a spatially 
periodic diffusion coefficient with small spatial period e, while the macroscopic (effective) 
equation describes the averaged behaviour. In the limit of e going to zero, this effective 
equation is the classical homogenized equation. Our goal is to approximate the effective 
equation by using only the microscopic equation in a set of small boxes. In [29], we already 
studied the gap-tooth scheme for periodic reaction-diffusion homogenization problems. We 
showed that the gap-tooth scheme approximates a finite difference scheme for the homoge- 
nized equation, when the averaged gradient is constrained at the box boundaries. However, 
generally, a given microscopic code only allows us to run with a set of predefined boundary 
conditions. It is highly non-trivial to impose macroscopically inspired boundary conditions 
on such microscopic codes, see e.g. [23] for a control-based strategy. Therefore, we cir- 
cumvent this problem here by introducing buffer regions at the boundary of each small 
box, which shield the short-term dynamics within the computational domain of interest 
from boundary effects. One then uses the microscopic code with its built-in boundary con- 
ditions. In this paper, we study the resulting gap-tooth scheme with buffers, which was 
introduced in [28, 29], when used inside a patch dynamics scheme, and analyze the relation 
between buffer size, time step and accuracy for a model diffusion homogenization problem. 
The analysis in this context is important, because we can clearly show the influence of 
the microscopic scales on the accuracy of the solution for this model problem. However, 
we emphasize that the real advantage of the method lies in its applicability for non-PDE 
microscopic simulators, e.g. kinetic Monte Carlo or molecular dynamics. 

It is worth mentioning that many numerical schemes have been devised for the ho- 
mogenization problem. Hou and Wu developed the multi-scale finite element method that 
uses special basis functions to capture the correct microscopic behaviour [13, 14]. Schwab, 
Matache and Babuska have devised a generalized FEM method based on a two-scale finite 
element space [26, 30]. Runborg et al. [27] proposed a time-stepper based method that 
obtains the effective behaviour through short bursts of detailed simulations appropriately 
averaged over many shifted initial conditions. The simulations were performed over the 
whole domain, but the notion of effective behaviour is identical. In their recent work, E and 
Engquist and collaborators address the same problem of simulating only the macroscopic 
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behaviour of a multiscale model, see e.g. [1, 7]. In what they call the heterogeneous mul- 
tiscale method, a macroscale solver is combined with an estimator for quantities that are 
unknown because the macroscopic equation is not available. This estimator subsequently 
uses appropriately constrained runs of the microscopic model [7]. It should be clear that 
patch dynamics does exactly this: by taking a few gap-tooth steps, we estimate the time 
derivative of the unknown effective equation, and give this as input to an ODE solver, 
such as projective integration. The difference in their work is that, for conservation laws, 
the macro-field time derivative is estimated from the flux of the conserved quantity; their 
generalized Godunov scheme is based on this principle. 

This paper is organized as follows. In section 2, we describe the model homogenization 
problem. In section 3, we show how to use the gap-tooth scheme to approximate the time 
derivative of the unavailable macroscopic equation. We prove a consistency result and 
propose a simple heuristic to obtain a sufficient buffer size. We also discuss to what extent 
the results depend on the specific setting of our model problem. In section 4, we describe the 
full patch dynamics algorithm and give some comments on stability. Section 5 contains some 
numerical examples which illustrate the accuracy and efficiency of the proposed method, 
and we conclude in section 6. 

2 The homogenization problem 

As a model problem, we consider the following parabolic partial differential equation, 



where a(y) = a (x/e) is uniformly elliptic and periodic in y and e is a small parameter. We 
choose homogeneous Dirichlet boundary conditions for simplicity. 

According to classical homogenization theory [5], the solution to (1) can be written as 
an asymptotic expansion in e, 



where the functions Ui(x, y, t) = Ui(x,x/e, t), i = 1, 2, . . . are periodic in y. Here, uo(x, t) is 
the solution of the homogenized equation 



d t u e (x, t) = d x (a (x/e) d x u € (x, t)) , in [0, T) x [0, 1] 
u e (x,0) = u°(x) eL 2 ([0, 1]), u € (0,t) = u e (l,t) =0, 



(1) 



oo 




(2) 



1=1 



d t u (x,t) = d x (a*d x u (x,t)) , in [0,T) x [0, 1] 
u {x, 0) = u°(x) G L 2 ([0, 1]), u (0, t) = «o(l, t) = 0, 



(3) 



the coefficient a* is the constant effective coefficient, given by 




(4) 
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and x(y) is the periodic solution of 

ly ( a(y) o| X(y) ) = a^ a(y) ' (5) 

the so-called cell problem. The solution of (5) is only defined up to an additive constant, so 
we impose the extra condition 

1 

X{y)dy = 0. 



o 



From this cell problem, we can derive m(x,y,t) = d x uox(y)- We note that in one space 
dimension, an explicit formula is known for a*, [5], 



1 1 



a (y) 



dy 



-i 



(6) 



These asymptotic expansions have been rigorously justified in the classical book [5], see 
also [6]. Under the assumptions made on a(x/e), one obtains strong convergence of u t (x,t) 
to uo(x,t) as e — > in L 2 ([0, 1]) x C([0,T)). Indeed, we can write 

\\u e (x,t) - u (x,t)\\ L2{[0yl]) < C e, (7) 

uniformly in t. 

It is important to note that the gradient of u(x, t) is given by 

d x u t (x,t) = d x u (x,t) + d y ux(x,y,t) + O(e), (8) 

from which it is clear that the micro-scale fluctuations have a strong effect on the local 
detailed gradient. 

Using the gap-tooth scheme, we will approximate the homogenized solution uo(x,t) by 
a local spatial average, defined as 

^ px+h/2 

U(x,t) = S h (u € ){x,t) = - u e {£,t)d£. 

n Jx-h/2 

It can easily be seen that that U(x,t) is a good approximation to uo(x,t) in the following 
sense. 

Lemma 2.1. Consider u e (x,t) to be the solution of (1), and uo(x,t) to be the solution of 
the associated homogenized equation (3). Then, assuming 

h = 0(eP), p 6(0,1), (9) 

the difference between the homogenized solution uo(x,t) and the averaged solution U(x,t) is 
bounded by 

\\U(x,t)-u (x,t)\\ Laam]) <dh 2 + C 2 e 1 -v. (10) 

For a proof, we refer to [29, Lemma 3.1]. Note that this error bound can be improved 
if we have more knowledge about the convergence of u e to uq (e.g. in L°°([0, 1])). 
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3 Estimation of the time derivative 



We devise a scheme for the evolution of the averaged behaviour U(x,t), while making only 
use of the given detailed equation (1). Moreover, we assume that a time integration code 
for (1) has already been written and is available with a number of standard boundary 
conditions, such as no- flux or Dirichlet. We also assume that the order d of the unavailable 
macroscopic equation (the highest spatial derivative) is known. A strategy to obtain this 
information is given in [22]. So, we know that the macroscopic equation is of the form 

d t u = F(u,d x u,...,d%u,t), (ii) 

where dt denotes the time derivative and d% denotes the fe-th spatial derivative. 
3.1 The gap-tooth scheme with buffers 

Suppose we want to obtain the solution of (11) on the interval [0, 1], using an equidistant, 
macroscopic mesh Tl(Ax) := {0 = xo < x\ = xq + Ax < . . . < xn = 1}- For convenience, 
we define a macroscopic comparison scheme, which is a space-time discretization for (11) in 
the assumption that this equation is known. We will denote the numerical solution of this 
scheme by ~ U(xi,t n ). Here, we choose as a comparison scheme a forward Euler/spatial 
finite difference scheme, which is defined by 

jjn+6 = S (U n ,t n ; St) = U n + 5t F{U n , D^U* 1 ), . . . , D d {U n ), t n ), (12) 

where D k (U n ) denotes a suitable finite difference approximation for the k-th spatial deriva- 
tive. 

Since equation (11) is not known explicitly, we construct a gap-tooth scheme to approx- 
imate the comparison scheme (12). We denote the solution of the gap-tooth scheme by 
UT 1 ss [/"". The gap-tooth scheme is now constructed as follows. Consider a small interval 
(box, tooth) of length h around each mesh point, as well as a larger buffer interval of size 
H > h. (See figure 1.) We will perform a time integration using the microscopic model 
(1) in each box of size H, and we provide this simulation with the following initial and 
boundary conditions. 

Initial condition We define the initial condition by constructing a local Taylor expansion, 
based on the (given) box averages U™, i = 0, . . . , N, at mesh point Xi and time t n , 

u*(x,t n ) = Y J DHu n ) [X ~ k f\ x £ [xi — *, Xi + |], (13) 
fc=o 

where d is the order of the macroscopic equation. The coefficients D^(U n ), k > are the 
same finite difference approximations for the k-th spatial derivative that would be used in 
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Figure 1: A schematic representation of the gap-tooth scheme with buffer boxes. We choose a number 
of boxes of size h around each macroscopic mesh point xi and define a local Taylor approximation 
as initial condition in each box. Simulation is performed inside the larger (buffer) boxes of size H, 
where some boundary conditions are imposed. 



the comparison scheme (12), whereas D®(U n ) is chosen such that 

^ i-Xi+h/2 



u*(£,i„)d£ = tff. (14) 

n Jxi-h/2 

For example, when d = 2, and using standard (second-order) central differences, we have 
fj n — r >fj n _i_ fj n fj n — fj n u 2 

DW n ) = 1+1 £ 2 Dht/n) = Ax = ^-Y2 mn) - (15) 

The resulting initial condition was used in [29], where it was derived as an interpolating 
polynomial for the box averages. 



Boundary conditions The time integration of the microscopic model in each box should 
provide information on the evolution of the global problem at that location in space. It is 
therefore crucial that the boundary conditions are chosen such that the solution inside each 
box evolves as if it were embedded in the larger domain. We already mentioned that, in many 
cases, it is not possible or convenient to impose macroscopically-inspired constraints on the 
microscopic model (e.g. as boundary conditions). However, we can introduce a larger box of 
size H > h around each macroscopic mesh point, but still only use (for macro-purposes) the 
evolution over the smaller, inner box. The simulation can subsequently be performed using 
any of the built-in boundary conditions of the microscopic code. Lifting and (short-term) 
evolution (using arbitrary available boundary conditions) are performed in the larger box; 
yet the restriction is done by processing the solution (here taking its average) over the inner, 
small box only. The goal of the additional computational domains, the buffers, is to buffer 
the solution inside the small box from the artificial disturbance caused by the (repeatedly 
updated) boundary conditions. This can be accomplished over short enough time intervals, 
provided the buffers are large enough] analyzing the method is tantamount to making these 
statements quantitative. 
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The idea of a buffer region was also introduced in the multiscale finite element method of 
Hou (oversampling) [13] to eliminate boundary layer effects; also Hadjiconstantinou makes 
use of overlap regions to couple a particle method with a continuum code [12]. If the 
microscopic code allows a choice of different types of microscopic boundary conditions, 
selecting the size of the buffer may also depend on this choice. 



The algorithm The complete gap-tooth algorithm to proceed from t n to t n + St is given 
below: 

1. Lifting At time t n , construct the initial condition u l (x,t n ), i = 0, ... ,N using the 
box averages U™, as defined in (13). 

2. Simulation Compute the box solution u l {x,t), t > t n , by solving equation (1) in the 
interval [ajj — H/2, Xi + H/2] with some boundary conditions up to time t n +s = t n + 5t. 
The boundary conditions can be anything that the microscopic code allows. 

3. Restriction Compute the average C/ J n+<5 = 1/h J^-h/^ ^*(£> *n+5)d£ over the inner, 
small box only. 

It is clear that this procedure amounts to a map of the macroscopic variables U n at 
time t n to the macroscopic variables at time t n +5, i-e. a "coarse to coarse" time 5i-map. We 
write this map as follows, 

fjn+s = s d (u n : tn ; St, H) = U n + St F d (U n , t n ; St, H), (16) 
where we introduced the time derivative estimator 

fjn+5 _ fjn 

F d (U n ,t n ;St,H) = . (17) 

The superscript d denotes the highest spatial derivative that has been prescribed by the 
initialization scheme (13). The accuracy of this estimate depends on the buffer size H, the 
box size h and the time step St. 



3.2 Consistency 

To analyze convergence, we solve the detailed problem approximately in each box. Because 
/i > e, we can resort to the homogenized solution, and bound the error using equation (7). 
It is important to note that we use the homogenized equation for analysis purposes only. 
The algorithm uses box averages of solutions of the detailed problem (1), so it does not 
exploit more knowledge of the homogenized equation than the order d. We choose to study 
convergence in the case of Dirichlet boundary conditions, but we will show numerically that 
the results do not depend crucially on the type of boundary conditions. 
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We first relate the gap-tooth time-stepper as constructed in section 3.1 to a gap-tooth 
time-stepper in which the microscopic equation has been replaced by the homogenized 
equation. 

Lemma 3.1. Consider the model equation, 

d t u e (x,t) = d x (a (x/e) d x u e (x,t)) , (18) 

where a(y) = a (x/e) is periodic in y and e <C 1, with initial condition u e (x,0) = u (x) and 
Dirichlet boundary conditions 

u € {-H/2, t) = u°(-H/2), u e (H/2, t) = u°(H/2). (19) 
For e — * 0, this problem converges to the homogenized problem 

d t u (x, t) = d x (a*d x u (x, t)) (20) 
with initial condition Uq(x, 0) = u°(x) and Dirichlet boundary conditions 

u € {-H/2, t) = u°(-H/2), u e (H/2, t) = u°(H/2). (21) 

and the solution of (18)-(19) converges pointwise to the solution of (20)- (21), with the 
following error estimate 



\\u e (x,t) - u (x,t)\\ L 2 i{ _ H/2tH/2]) < C 3 e. 

This is a standard result, whose proof can be found in e.g. [2, 6]. 
We now define two gap-tooth time-steppers. Let 

U n+S = S 2 (U n , t n - St, H) = U n + St F 2 (U n , t n ; St, H) 



(22) 



(23) 



be a gap-tooth time-stepper that uses the detailed, homogenization problem (18)-(19) inside 
each box, and 



U 



n+5 



S 2 (U n , t n ; St, H) = U n + St F\U n , t n ; St, H) 



(24) 



be a gap-tooth time-stepper where the homogenization problem for each box has been 
replaced by the homogenized equation (20)-(21). The box initialization is done using a 
quadratic polynomial as defined in (15). 

We can apply [29, Lemma 4.2] to bound the difference between F 2 (U ,t n ; St, H) and 
F 2 (U,t n ;5t,H). 

Lemma 3.2. Consider U n+S = S 2 (U n ,t n ;St,H) and U n+S = S 2 {U n ,t n ;St,H) as defined 
in (23) and (24), respectively. Assuming U n = U n , h = 0(e p ), p G (0, 1), e — ► 0, we have 



01 



U! 



n+5 



< C^-vl 2 , 



and therefore 



F 2 {U n ,t n ; St, H) - F 2 (U n ,t n ; St, H] 



< Ca 



E l-p/2 

St ' 
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Again, note that the error estimate can be made sharper if additional knowledge of the 
convergence of u e to no is available. 

It can easily be checked that the averaged solution U(x, t) also satisfies the diffusion 
equation (20). Therefore, we define the comparison scheme (12) for the model problem as 

jjn+s = S (U n ,t n ;5t) 

= U n + 5t F{U n ,D 1 {U n ),D 2 {U n ),t n ) 

= U n + St[a* D 2 {U n )]. (25) 

The following theorem compares the gap-tooth time derivative estimator F 2 (U n ,t n ; 5t, H) 
with the finite difference time derivative used in (25). 

Theorem 3.3. Consider the gap-tooth time- stepper for the homogenized equation, as defined 
by (24-), and the corresponding comparison scheme (25). Assuming U n = U n , and defining 
the error 

E(5t,H)= F 2 (U n ,t n ;5t,H) -a* D 2 {U n ) 
we have the following result for 5t/H 2 — *■ 0, h <C H , 

h 2 \ A , * 2 St 



E(5t, H) < Id + C 2 - \l - exp(-aV-p) \ (26) 

Proof. First, we solve the equation (20)-(21) analytically inside each box, with initial con- 
dition given by (13). Using the technique of separation of variables, we obtain 

u\x, t) =&r - ^D 2 (U n ) + D 2 n )^- + D\(U n ){x - Xi ) 

E, / ^rn 2 TT 2 A . f m-K . H s 

a m expl -a ——{t-tn) J sm I — {x - x t - — J 

where 

2 f x >+ H / 2 1 „ * / H 2 \ (m-K, . 

m = Hj x _ H2 2 Di{u ) \ {x ~ Xi) ~ t) sin [ir ix " Xi ~ ~ ] ' 

This can be simplified to 

2H 2 D 2 (U n ) ((-1)™ - 1) 



m 3 7r 3 



which yields the following solution, 
u\x, t) = U?- ^D 2 (U n ) + D 2 (U n )^- + Dj(U n )(x - Xl ) 



m=l v 7 v 7 
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When taking the average over a box of size ft, we obtain, 

1 /• '• • - ( . ^ ft 2 „o,^, „o,^,H 2 

'xi-h/2 



/ u\x,t)dx = U? - Df(U n ) + Df(U n ) — 



m 

with a m determined by 



1 r^+V2 /(2m-l)7r / H. . 

a m = - I sin — (x - Xj - — J dx 



ft Jxi-h/2 \ H 2 

' cos — — (JT + ft) - cos - ft) 



(2m - l)ftvr V V 2H K 'J \ 2H 

V ; (2m - l)ft7r V 2 # 
The coefficients a m tend to 1 in absolute value as ft — > 0. To obtain the time derivative 
estimate F(U n ,t n ;5t,H), we proceed as follows: 

i rXi+h/2 

F?(U n ,t n ;6t,H) = — u\x,t n + 8t)-u l {x,t n )dx 

ot n J Xi -h/2 

1^ 4jff 2 / ,(2m-l) 2 vr 2 \ \ 

= £ E ( 2m - 1)3,3 A 2 (^) (exp ^ g2 J St) ~ lj 

m=l x x / ' 

AB^tW- ( Um (exp(-^(2m-l) 2 vr 2 g)-l) \ 
H^i (2m-l)3vr3 J 

where we introduced £ = 5t/H 2 . It can easily be checked (e.g. using Maple) that 

\hnF?(U n ,t n ;5t,H) = a*Df(U n ), 

which already shows that the gap-tooth scheme is consistent in this limit. Obtaining an 
error bound in terms of £ is somewhat more involved. We split F?(U n ,t n ; 5t, H) as follows, 

F(U n ,t n ;6t,H)=F 1 +F 2 , 

with F\ and F 2 defined as 

F AT*(T>\ 1 ( V (~l) m (exp(-a*(2m-l) 2 ^)-l) \ 

p 2 = 4^(^)1 f v (^-(-m(«p(-^-iM)-i)) . 

£ \^ (2m - 1)3.3 y 

We now show that F\ approaches the correct estimate exponentially. Some algebraic ma- 
nipulation results in 

A-^(-) « ^ (g ( - ir (6XP fcf V } ^ " ^ ) - ^ 



4n2 ™f,_,™ + i / 1 - a*(2m - l) 2 vr 2 £ - exp (-q'(2m - l) 2 vr 2 g) 
^ V (2m-l)3vr3£ 
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Therefore, we have 



A -a*Df(U n )\\ = 

AT*(T>\ ff_,F+i„« A-a*(2m- l) 2 ^ 2 g ~ exp (-g*(2m - l) 2 ^) 



m=l 



< c 



'(2m - 1) 3 tt 3 £ 
1 - a*7r 2 £ - exp (-a*7r 2 ^) 



a*7r 2 £ 



< C(l-exp(-aV£)) (29) 
It remains to show the asymptotic behaviour of F2. 



4A 2 (^)E(- 1 )" A ^ ^ fexp (-a* (2m - l) 2 ^) - 1 



m=l 



(2m - 1) 3 tt 3 



< C rMMih * j (1 - exp (-aV£)) 

fr 2 1 - exp(-a*7r 2 g) 
" H 2 £vr 3 

< C ^(l-exp(-aVC)) (30) 
The combination of (29) and (30) proves the theorem. □ 

The error bound (26) clearly shows an exponential decay of the error as a function of 
5t/H 2 when the microscopic problem is replaced by the effective equation. The restriction 
(here performed by taking the box average) also affects the accuracy of the estimate. Ideally, 
one would just use the effective function value at x = X{ inside each box (this corresponds 
to h = 0), but when microscopic scales are present, this value is generally impossible to 
obtain. 

We illustrate this result numerically. 

Example 3.4. Consider the model problem (20) with a* = 0.45825686 as a microscopic 
problem on the domain [0, 1] with homogeneous Dirichlet boundary conditions and initial 
condition u(x, 0) = 1 — 4(x — 1/2) 2 . To solve this microscopic problem, we use a second order 
finite difference discretization with mesh width Sx = 2 • 10 -7 and lsode as time-stepper. 
The concrete gap-tooth scheme for this example is defined by the initialization (15). We 
compare a gap-tooth step with h = 2 ■ 10~ 3 and Ax = 1 • 10 _1 with the reference estimator 
a*D 2 (U n ). Figure 2 shows the error with respect to the finite difference time derivative as a 
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Figure 2: Error of the gap-tooth estimator F 2 (U n ,t n ; 5t, H) (which uses the homogenized problem 
(20)-(21) inside each box) with respect to the finite difference time derivative a*D 2 (U n ) on the same 
mesh. Left: Error with respect to H for fixed St. Right: Error with respect to St for fixed H. 

function of H (left) and St (right). It is clear the convergence is in agreement with theorem 
3.3. The stagnation for large buffer sizes is due to the finite accuracy of the microscopic 
solver. 

We are now ready to state the general consistency result. 

Theorem 3.5. Let U n+S = S 2 (U n ,t n ;S t , H) be a gap-tooth time-stepper for the homoge- 
nization problem (18)-(19), as defined in (23), and U n+S = S(U n , t n ; 5t) a comparison finite 
difference scheme as defined in (25). Then, assuming U n = U n , we have, 

f X ~vl 2 / h 2 \ / St 

\\F 2 (U n ,t n ;5t,H)-a*D 2 (U n )\\<C 4 ^— +C 5 1 + — 1 - exp(-aV-^) 



St J \~ ~"^ v ~ " H 2 

v " V v 

averaging boundary conditions 

(31) 

Proof. This simply follows by combining theorem 3.3 with lemma 3.1. □ 

Formula (31) shows the main consistency properties of the gap-tooth estimator. The 
error decays exponentially as a function of buffer size, but the optimal accuracy of the 
estimator is limited by the presence of the microscopic scales. Therefore, we need to make 
a trade-off to determine an optimal choice for H and 5t. The smaller St, the smaller H can 
be to reach optimal accuracy (and thus the smaller the compational cost), but smaller St 
implies a larger optimal error. This is illustrated in the following numerical example. 

Example 3.6. Consider the model problem (18) with 

a(x/e) = 1.1 + sin(2vrx/e), e = 1 • 10~ 5 (32) 
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Figure 3: Error of the gap-tooth estimator F(U n , t n ; St, H) (which uses the detailed, homogenization 
problem (18)-(19) inside each box) with respect to the finite difference time derivative a*D 2 (U n ) 
on the same mesh. Left: Error with respect to H for fixed St. Right: Error with respect to St with 
fixed H. 



as a microscopic problem on the domain [0, 1] with homogeneous Dirichlet boundary condi- 
tions and initial condition u(x, 0) = 1 — 4(x — 1/2) 2 . This diffusion coefficient has also been 
used as a model example in [1, 29]. To solve this microscopic problem, we use a second order 
finite difference discretization with mesh width 8x = 1 • 10~ 7 and lsode as time-stepper. 
The concrete gap-tooth scheme for this example is defined by the initialization (15). We 
compare a gap-tooth step with h = 2 • 10~ 3 and Ax = 1 • 10" 1 with the reference estima- 
tor a*D 2 (U n ), in which the effective diffusion coefficient is known to be a* = 0.45825686. 
Figure 3 shows the error with respect to the finite difference time derivative as a function 
of H (left) and 5t (right). It is clear that the convergence is in agreement with theorem 
3.5. We see that smaller values of 5t result in larger values for the optimal error, but the 
convergence towards this optimal error is faster. 



3.3 Choosing the method parameters 

When performing time integration using patch dynamics, one must determine a macroscopic 
mesh width Ax, an inner box size h, a buffer box size H and a time step 5t. These method 
parameters need to be chosen adequately to ensure an accurate result. Since the gap-tooth 
estimator approximates the time derivative that would be obtained through a method- 
of-lines discretization of the macroscopic equation, the macroscopic mesh width Ax can 
be determined by macroscopic properties of the solution only, enabling reuse of existing 
remeshing techniques for PDEs. The box width h has to be sufficiently large to capture 
all small scale effects, but small enough to ensure a good spatial resolution. Here, we just 
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choose h S> e. In our simplified setting, where the microscopic model is also a partial 
differential equation, we are free to choose 8t, which allows us to illustrate the convergence 
properties of the method. However, in practical problems, the choice of 5t will be problem- 
dependent, since it will need to be chosen large enough to deduce reliable information on 
the macroscopic time derivative. 

Therefore, we focus on determining the buffer width H, assuming that all other pa- 
rameters have already been fixed. From theorem 3.5, it follows that the desired value of 
H depends on the effective diffusion coefficient a*, which is unknown. We thus need to 
resort to a heuristic. Consider the model problem (18)-(19) inside one box, centered around 
x = 5-10" 1 , with ff = 8-10" 3 , with initial condition u°(x) = 1 -4(x- 1/2) 2 . The diffusion 
coefficient is given by (32), see example 3.6. Denote the solution of this problem by u(x,t), 
and define 



F(x,t) = S h (u{x,t) -u(x,0)) 

1 ft- X0+h / 2 u(£,t)-u(Z,0) 



t ■lt-xo-h/2 h 



d£, (33) 



with h = 2 ■ 1CT 3 and x € [(-H + h)/2, (H - h)/2]. Figure 4 (left) shows F(x,t) for a 
number of values of t. We clearly see how the error in the estimator propagates inwards 
from the boundaries. The same function is plotted on the right, only now the microscopic 
model is the reaction-diffusion equation 

d t u e (x, t) = d x (a(x/e)d x u e (x, t)) + u e (x, t) ( 1 



1.2 + sin(27nr) / (34) 
u e (-H/2,t) = u°(-H/2), u e (H/2,t) =u°{H/2), 

again with a(x/e) defined as in (32). In the presence of reaction terms, F(x,t) is no longer 
constant in the internal region. Based on these observations, we propose the following test 
for the quality of the buffer size, 

\\F(0,5t) - F(0,(1 - a)5t)\\ <Tz, < a < 1. (35) 

Figure 5 shows this heuristic, together with the error, as a function of H for 5t = 5 • 10~ 6 
and a = 0.04. It is clear that the computed quantity in (35) is proportional to the error for 
sufficiently large H. However, this heuristic is far from perfect, since the simulations inside 
each box can converge to a steady state due to the Dirichlet boundary conditions. If this 
steady state is reached in a time interval smaller than St, equation (35) will underestimate 
the error, resulting in an insufficient buffer size H getting accepted. However, as soon as 
the problem-dependent parameters a and Tr have been determined, this heuristic can be 
used during the simulation to check whether the currently used buffer size is still sufficient. 



15 



-10 1 1 ' • ' 1 1 1 1 0.7 1 ■ ' • ' 1 1 1 1 

-0.004 -0.002 0.002 0.004 -0.004 -0.002 0.002 0.004 

X - ,T X - ,T 

Figure 4: The function F(x, t) as denned in equation (33) for a number of values of time, using a 
buffer size H = 8 • 1CT 3 and h = 2 • 1CT 3 . Left: the model diffusion problem (18-19). Right: the 
reaction-diffusion equation 34. The estimate clearly gets affected by the boundary conditions as 
time advances. 
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Figure 5: Error of the gap-tooth estimator (dashed) and heuristic error estimate (solid) as a function 
of buffer size for the model equation (18) with diffusion coefficient (32) for St = 5 ■ 10~ 6 and a = 0.04. 
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Figure 6: Error of the gap-tooth estimator F(U n , t n ; St, H) (using the microscopic problem (36) with 
diffusion coefficient (32) in each box) with respect to the finite difference time derivative a*D 2 (U n ) 
on the same mesh. 



3.4 Discussion 



Other boundary conditions In section 3.2, we studied the convergence of the gap-tooth 
estimator both analytically and numerically in the case of Dirichlet boundary conditions. 
We will now show numerically that the results obtained in that section do not depend 
crucially on the type of boundary conditions. Consider again the diffusion problem (18), 
with the diffusion coefficient defined as in (32), see also example 3.6. We construct the gap- 
tooth time derivative estimator F(U n ,t n ;5t,H) as outlined in section 3.1, but now we use 
no-flux instead of Dirichlet boundary conditions. In each box, we then solve the following 
problem, 

d t u t (x,t) = d x (a(x/e)d x u e (x,t)) , 

(36) 

d x u e (-H/2, t) = 0, d x u t (H/2, t)=0 

The concrete gap-tooth scheme that is used, as well as the corresponding finite difference 
comparison scheme, are defined by the initialization (15). Figure 6 shows the error with 
respect to the finite difference time derivative a*D 2 (U n ). We see qualitatively the same 
behaviour as for Dirichlet boundary conditions. 

In general, the choice of boundary conditions might influence the required buffer size. In 
the ideal case, where the boundary conditions are chosen to correctly mimic the behaviour 
in the full domain, we can choose H = h. Then there is no buffer and the computational 
complexity is, in some sense, optimal. For reaction-diffusion homogenization problems, 
this can be achieved by constraining the averaged gradient around each box edge [29]. In 
situations where the correct boundary conditions are not known, or prove impossible to 
implement, one is forced to resort to the use of buffers. 
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Microscopic simulators It is possible that the microscopic model is not a partial dif- 
ferential equation, but some microscopic simulator, e.g. kinetic Monte Carlo or molecular 
dynamics code. In fact, this is the case where we expect our method to be most useful. 
In this case, the lifting step, i.e. the construction of box initial conditions, becomes more 
involved. In general, the microscopic model will have many more degrees of freedom, the 
higher order moments of the evolving distribution. These will quickly become slaved to the 
governing moments (the ones where the lifting is conditioned upon), see e.g. [19, 24]. The 
crucial assumption in theorem 3.5 is that the solution in each box evolves according to the 
macroscopic equation. For a microscopic simulation, this will usually mean that we need to 
construct an initial condition in which, for example, a number of higher order moments are 
already slaved to the governing moments (so-called mature initial conditions). To this end, 
it is possible to perform a constrained simulation before initialization to create such mature 
initial conditions [9, 15]. If this is not done, the resulting evolution may be far from what 
is expected, see [21] for an illustration in the case of a lattice-Boltzmann model. 



4 Patch dynamics 

Once a good gap-tooth time derivative estimator has been constructed, it can be used as 
a method-of-lines spatial discretization in conjunction with any time integration scheme. 
Consider for concreteness the forward Euler scheme for (11), given by 

jjn+i = r/« + AfF(C/ n ,D 1 (C/ n ) ) ...,D d (C/ 7l ),t n ), (37) 
which we will abbreviate as 

jjn+l =U n + At F(U n ,t n ) (38) 

and the corresponding patch dynamics scheme 

U n+1 = U n + At F d (U n ,t n ;5t,H), (39) 

where F d (U n , t n ; St, H) is defined as in (17). Theorem 3.5 establishes the consistency of the 
gap-tooth estimator. In order to obtain convergence, we also need to prove stability. For 
this purpose, we define the class K of discrete functions with bounded divided differences 
up to order d on the numerical grid (xj, t n ), i = 0, . . . , JV; t n = nAt, n = 0, . . . , T/ At, as 

K = {{U?}\\\D% X U? || < C a for a < d, nAt < T} , 

where d is the highest spatial derivative present in equation (11), D% x is the finite difference 
operator of order a on a mesh of width Ax, and C a is independent of At and Ax. We can 
then make use of [7, Theorem 5.5] to state the following result. (See also [29, Theorem 4.5] 
in the case of constrained gradient boundary conditions.) 
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Theorem 4.1. Consider the patch dynamics scheme (39) and the corresponding finite dif- 
ference comparison scheme (37). Assume that {U n }, {U n }, {U n } G K, U° = U° = U° and 
the comparison scheme (37) is stable, then we have 

\\U n -u (xi,t n )\\ < d(Ax k + At) + C 2 max \\F d (U n , t n ; St, H) - F(U n , t n )\\. 

0<k<T/At 

in the limit where 

\\F d (U n ,t n ;5t,H)-F(U n ,t n )\\ 0. 

Thus the patch dynamics scheme is stable if the finite difference comparison scheme is 
stable. We note that, although this result is very general, its applicability is limited due to 
the assumption that U n E K, which has to be checked separately. Therefore, this result does 
not prevent the patch dynamics scheme to become unstable, due to e.g. an insufficient buffer 
size H. However, we can study the stability of the patch dynamics scheme numerically by 
computing the eigenvalues of the time derivative estimator as a function of H. 

Consider the homogenization diffusion equation (18) with the diffusion coefficient a(x/e) 
given by (32). The homogenized equation is given by (20) with a* = 0.45825686. In this 
case, the time derivative operator F(U n ,t n ) in the comparison scheme (37) has eigenvalues 

4a* 

X k = ~ A-2 ^ k AX ) ' ( 40 ) 

Ax z 

which, using the forward Euler scheme as time-stepper, results in the stability condition 

At 1 

maxll + AfcAtt < 1 or — — ~ < -a* 
k Ax z 2 

It can easily be checked that the operator F(U n ,t n ; 5t, H) is linear, so we can interpret the 
evaluation of F(U n , t n ; 5t, H) as a matrix-vector product. We can therefore use any matrix- 
free linear algebra technique to compute the eigenvalues of F(U n ,t n ;5t,H), e.g. Arnoldi. 
We choose to compute F(U n ,t n ;St,H) and F(U n ,t n ) on the domain [0,1] with Dirichlet 
boundary conditions, on a mesh of width Ax = 0.05 and with an inner box width of 
h = 2 • 10 -3 . We choose 8t = 5 • 10 -6 and compute the eigenvalues of F(U n ,t n ;5t, H) as 
a function of H. The results are shown in figure 7. Two conclusions are apparent: since 
the most negative eigenvalue for F(U n ,t n ;5t,H) is always smaller in absolute value than 
the corresponding eigenvalue of F(U n ,t n ) the patch dynamics scheme is always stable if 
the comparison scheme is stable. Moreover, we see that, with increasing buffer size H, the 
eigenvalues of F(U n ,t n ; St, H) approximate those of F(U n ,t n ), which is an indication of 
consistency. 



5 Numerical results 

We will consider two example systems to illustrate the method. The first example is a system 
of two coupled reaction-diffusion equations, which models CO oxidation on a heterogeneous 
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Figure 7: Spectrum of the estimator F(U n ,t n ;5t,H) (dashed) for the model equation (18) with 
diffusion coefficient (32) for H = 2 ■ 10 -3 , 4 • 10 -3 , . . . , 2 • 10 -2 and 5t = 5- 10~ 6 , and the eigenvalues 
(40) of F(U n ,t n ) (solid). 



catalytic surface. Due to the reaction term, the proof of theorem 3.5 is strictly speaking not 
valid, but nevertheless the conclusions are the same. The second example is the Kuramoto- 
Sivashinsky equation. This fourth-order non-linear parabolic equation is widely used e.g. in 
combustion modeling. The patch dynamics scheme with buffers also works in this case, 
showing the more general applicability of the method. All computations were performed in 
Python, making use of the SciPy package [16] for scientific computing. 



5.1 Example 1: A nonlinear travelling wave in a heterogeneous excitable 
medium 



Consider the following system of two coupled reaction-diffusion equations, 

w(x, t) + b(x) 
a(x) 



dtu(x, t) = d%.u(x, t) + -u(x, t)(l — u(x, t)) ( u(x, t) 
d t w(x, t) = g(u(x, t)) - w(x, t), 



with 



0, u < 1/3, 
l-6.75u(l-u) 2 , 1/3<u<1, 

1, u>l. 



(41) 



(42) 



This equation models the spatiotemporal dynamics of CO oxidation on microstructured cat- 
alysts, which consist of, say, alternating stripes of two different catalysts, such as platinum, 
Pt, and palladium, Pd, or platinum and rhodium, Rh [11, 4, 32]. The goal is to improve 
the average reactivity or selectivity by combining the catalytic activities of the different 
metals, which are coupled through surface diffusion. In the above model, u corresponds to 
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the surface concentration of CO, w is a so-called surface reconstruction variable and g(u) 
is an experimentally fitted sigmoidal function. Details can be found in [17, 3]. 

In this model a and b and the time-scale ratio parameter 5 are physical parameters that 
incorporate the experimental conditions: partial pressures of O2 and CO in the gas phase, 
temperature, as well as kinetic constants for the surface. Here, we will study a domain of 
length L = 21 with a periodically varying medium: a striped surface that can be thought 
of as consisting of equal amounts of Pt and Rh, with stripe width e/2. The medium is then 
defined by 

a(x) = 0.84, b(x) = -0.025 + 0.725 sin(27rx/e), <5 = 0.025. (43) 

This particular choice of parameters is taken from [27], where an effective bifurcation anal- 
ysis for this model was presented. For these parameter values, the effective equation (given 
by (41)-(42) with 

a(x) = 0.84, b(x) = -0.025, 5 = 0.025, (44) 

supports travelling waves. It was shown in [27] that this conclusion remains true for the 
given heterogeneity. This was done by computing the effective behaviour as the average of a 
large number of spatially shifted realization of the wave. Here, using the gap-tooth scheme, 
the solution is spatially averaged inside each box, but the notion of effective behaviour is 
identical. We choose the small scale parameter e = 1 • 10 . 

The macroscopic comparison scheme for the effective equation (41-42)-(44) is defined as 
a standard second order central difference discretization in space on a macroscopic mesh of 
width Ax = 0.25, combined with a forward Euler time-stepper. The time-step is chosen as 
At = M0~ 2 , which ensures stability. The patch dynamics scheme for the detailed equation 
(41-43) is then obtained by using a gap-tooth estimator for the time derivative using the 
initialization (15) with the same forward Euler time-stepper. 



Accuracy We perform a numerical simulation for this model on the domain [0, L\ using 
the patch dynamics scheme. The gap-tooth parameters are given by h = 5 • 10~ 4 , H = 
1.5 • 10 -2 and St = 5 ■ 10 -7 . Inside each box, we used a finite difference approximation in 
space, with mesh width Sx = 1 ■ 10 -6 and lsode as time-stepper. The initial condition is 
given by 



u(x, 0) 



1, x € [8, 18] 
0, else 



w(x, 0) 



0.5 -0.05a?, 
0.07x - 0.46, 
-0.1a; + 2.6, 



x < 8, 

8 < x < 18, 

x > 18. 



The results are shown in figure 8. We clearly see both the initial transient and the final 
travelling wave solution. For comparison purposes, the same computation was performed 
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Figure 8: Left: solution of equation (41-43) using the patch dynamics scheme as a function of space 
and time. Colors indicate values (blue = 1, red = 0). Right: snapshots of the solution at certain 
moments in time, clearly showing the approach to a travelling wave solution. 



using the finite difference comparison scheme for the effective equation. We also computed 
an "exact" solution for the effective equation using a much finer grid (Ax = 5 • 10 -3 and 
At = 1-10~ 5 . Figure 9 shows the errors of the patch dynamics simulation with respect to the 
finite difference simulation of the effective equation and the "exact" solution, respectively. 
We clearly see that the patch dynamics scheme is a very good approximation of the finite 
difference scheme, and the error with respect to the exact solution is dominated by the error 
of the finite difference scheme. 

Efficiency Time integration using the patch dynamics scheme is more efficient than a 
complete simulation using the microscopic model, since the microscopic model is used only 
in small portions of the space-time domain (the patches). An obvious (but not always 
correct) way to study the efficiency is to compare the size of the total space-time domain 
with the size of the patches. In this example, the simulations are only performed in 6% 
of the spatial domain. Of course, when it is possible to apply physically correct boundary 
conditions around the inner box, the buffer boxes are not necessary, and the boxes would 
only cover 0.2 % of the space domain. For reaction-diffusion homogenization problems, we 
showed that buffer boxes are not required when we constrain the average gradient at the 
box boundary [29]. The gain in the spatial dimension is determined by the separation in 
spatial scales. It can be large when the macroscopic solution is smooth (few macroscopic 
mesh points are needed) and propagation of boundary artefacts is slow (small buffer box is 
sufficient). Note that in higher spatial dimensions, this gain can be even more spectacular. 

The gain in the temporal dimension can be determined similarly. In the example of 
section 5.1, the gap-tooth step was chosen as St = 5 • 10 -7 , whereas for macroscopic time 
integration, the forward Euler scheme was used with At = 1 • 1CP 2 . Therefore, in the 
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temporal dimension, we gain a factor of 2 • 10 . In more realistic applications, when the mi- 
croscopic model is not a partial differential equation, we expect this gain to be smaller, since 
additional computational effort will be required to remove the errors that were introduced 
during the lifting step, e.g. in the form of constrained simulation [9]. 

5.2 Example 2: Kuramoto— Sivashinsky equation 

Consider the Kuramoto-Sivashinsky equation 

dtu(x,t) = —ud x u(x,t) — d x u(x,t) — u(x,t)d x u(x,t), x G [0, 2tt], (45) 

with periodic boundary conditions. This equation is frequently used in the modelling of 
combustion and thin film flow. For the parameter value v = 4/15, it has been shown 
that the equation supports travelling wave solutions, see e.g. [20]. For the purpose of this 
example, both the microscopic and the macroscopic model are given by (45). 

To obtain the macroscopic comparison scheme, we discretize the second and fourth order 
spatial derivatives using second order central differences, on a macroscopic mesh of width 
Ax = 0.057T, combined with a forward Euler time integrator with time-step At = 1 • 10 -5 . 
This small macroscopic time-step arises due to the stiffness of the effective equation. We 
can accelerate time-stepping by wrapping a so-called projective integration method around 
the forward Euler scheme [8]. This scheme works as follows. First, we perform a number of 
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forward Euler steps, 

U k+1 ' N = U k > N + AtF(U k > N ,t n ), 
where, for consistency, U 0,N = U N , followed by a large extrapolation step 

U N+1 = (M + 1) U k+1 ' N - M U k ' N , M>k. 

Here, U k,N ~ U(N (M + k + 1) At + kAt). The parameters k and M determine the stability 
region of the resulting time-stepper. An analysis of these methods is given in [8]. It can be 
checked that, for this equation, choosing k = 2 and M = 7 results in a stable time-stepping 
scheme. 

The patch dynamics scheme is constructed by replacing the time derivative F(U k ' N ,t n ) 
by a gap-tooth estimator F 4 (U k,N ,t n ;St, H), obtained by the initialization (13), where 
we choose the order of the Taylor expansion to be d = 4. The coefficients D l k , k > are 
determined by the macroscopic comparison scheme. Inside each box, equation (45) is solved, 
on a mesh of width Sx = 1 • 10~ 5 , subject to Dirichlet and no-flux boundary conditions, 
using lsode as time-stepper. We fixed the box width h = 1 • 10 -3 . 

Consistency and efficiency Because of the fourth order term, theorem 3.5 is not 
proven. Therefore, we numerically check the consistency of the estimator, by computing 
the gap-tooth estimator F A (U k,N ,t n ;St, H) as a function of H for a range of values for St, 
and comparing the resulting estimate with the time derivative of the comparison scheme. 
As an initial condition, we choose u (x) = sin(2-7rx). The results are shown in figure 10 
(left). We see qualitatively the same behaviour as in section 3.2 for diffusion problems. 
There are two main differences. First, in this case the convergence is no longer monotonic, 
which explains the sharp peaks in the error curves. Also, because boundary artefacts travel 
inwards much faster, the gain will be much smaller. Indeed, the figure suggests that a 
good compromise between accuracy and efficiency would be to choose St = 4 • 10 -9 and 
H = 3tt ■ 10~ 2 . The reason for this behaviour is that the macroscopic equation contains 
reasonably fast time scales. Note that this is also the reason why the finite difference 
comparison scheme is forced to take small time-steps. 

For this choice of the parameters, the computations have to be performed in 60% of the 
spatial domain. However, for a forward Euler step, we only need to simulate in 1/25000 of 
the time-domain. Using the projective integration scheme therefore gives us a total gain 
factor of about 80000 in time. Again, we note that in real applications, this spectacular 
gain will partly be compensated by the additional computational effort that is required to 
create appropriate initial conditions. 

We can draw two main conclusions. The scheme allows to simulate higher order macro- 
scopic equations, and the gain in the space domain is heavily dependent on the separation 
of scales in the macroscopic equation. 
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Figure 10: Left: Error of the gap-tooth estimator F(U°'° ,t n ; 8t, H) with respect to the finite differ- 
ence time derivative F(U 0,0 ,t n ). Right: The function (33) as a function of x for a number of values 
of time. We clearly see how the estimate gets affected by the boundary conditions. 



Accuracy We perform a numerical simulation for this model on the domain [0, 2ir] using 
the patch dynamics scheme. The gap-tooth parameters are given by h = 1 • 10 -3 , H = 
3vr • 1CT 2 and 5t = 4 • 10" 9 . Inside each box, we used a finite difference approximation in 
space, with mesh width Sx = 1 • 10 -5 and lsode as time-stepper. The initial condition is 
given by 

-1, x E [0, 0.8-71"], 

u(x, 0) = <J -1 + 5(x - 0.8vr), x G [0.8vr, 1.5vr], 
2.5 - 7(x - 1.5vr), x G [i.5tt, 2vr], 

The results are shown in figure 11. We clearly see both the initial transient and the final 
travelling wave solution. For comparison purposes, the same computation was performed 
using the finite difference comparison scheme for the effective equation. Figure 12 shows 
the errors of the patch dynamics simulation with respect to the finite difference simulation. 
We see that during the transient phase the error oscillates somewhat, but once the trav- 
elling wave is steady the error increases linearly, due to a difference in the approximated 
propagation speed. Note that the error is significantly larger than for example 5.1, due to 
the fact that the estimator is less accurate, but also because the macroscopic time-step is 
much smaller, resulting in a larger number of estimations. 
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Figure 11: Left: solution of equation (45) using the patch dynamics scheme as a function of space 
and time. Colors indicate values (blue = 4, red = -4). Right: snapshots of the solution at certain 
moments in time, clearly showing the approach to a travelling wave solution. 
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Figure 12: Error of a patch dynamics simulation for equation (45) with respect to the a finite 
difference comparison scheme for the effective equation. We see that this error grows monotonic 
once the travelling wave has been reached, due to a slight difference in propagation speed. 
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6 Conclusions 



We described the patch dynamics scheme for multi-scale problems. This scheme approxi- 
mates an unavailable effective equation over macroscopic time and length scales, when only 
a microscopic evolution law is given; it only uses appropriately initialized simulations of 
the microscopic model over small subsets (patches) of the space-time domain. Because it is 
often not possible to impose macroscopically inspired boundary conditions on a microscopic 
simulation, we propose to use buffer regions around the patches, which temporarily shield 
the internal region of the patches from boundary artefacts. 

We analytically derived an error estimate for a model homogenization problem with 
Dirichlet boundary conditions. The numerical results show that the algorithm is more 
widely applicable. We showed the scheme is capable of giving good approximations for 
reaction-diffusion systems, as well as for fourth order PDEs, such as the Kuramoto-Sivashinsky 
equation. As such, these results are far more general than those of [29], which are restricted 
to reaction-diffusion problems due to the special choice of boundary conditions. 

We emphasize that, although analyzed for homogenization problems, the real advantage 
for the methods presented here lies in their applicability for microscopic models that are not 
PDEs, such as kinetic Monte Carlo, or molecular dynamics. Experiments in this direction 
are currently being pursued actively. 
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